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Abstract 

Many applications in genetic analyses utilize sampling distributions, which describe the 
probability of observing a sample of DNA sequences randomly drawn from a population. 
In the one-locus case with special models of mutation such as the infinite-alleles 
model or the finite-alleles parent-independent mutation model, closed-form sampling 
distributions under the coalescent have been known for many decades. However, no 
exact formula is currently known for more general models of mutation that are of 
biological interest. In this paper, models with finitely-many alleles are considered, 
and an van construction related to the coalescent is used to derive approximate closed- 
form sampling formulas for an arbitrary irreducible recurrent mutation model or for 
a reversible recurrent mutation model, depending on whether the number of distinct 
observed allele types is at most three or four, respectively. It is demonstrated empirically 
that the formulas derived here are highly accurate when the per-base mutation rate is 
low, which holds for many biological organisms. 
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1. Introduction 

An important problem in genetic analyses concerns computing the probability of observ- 
ing a randomly drawn sample of chromosomes under a given model of evolution. Popular 
applications of this probability computation include maximum likelihood estimation of model 
parameters and ancestral inference (see [19] for a nice introduction). The coalescent [14, 15] 
is a useful mathematical framework for performing model-based fuU-Ukelihood analyses, but 
in most cases it is intractable to obtain a closed-form formula for the probability of a given 
dataset. A well-known exception to this comphcation is the celebrated Ewens sampling for- 
mula (ESF) [3], which describes the stationary probability distribution of a sample configura- 
tion under the one-locus infinite-alleles model in the coalescent or the diffusion limit. A Polya- 
hke urn model interpretation [9] of the formula has been known for some time, and recently a 
new combinatorial proof of the ESF has been provided [6]. Furthermore, the ESF also arises 
in several interesting contexts outside biology, including random partition structures; the ESF 
is a special case of the two-parameter samphng formula [17, 18] for exchangeable random 
partitions. See [1] for examples of other interesting combinatorial connections. 
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In the case of finitely-many alleles, a closed-form sampling formula is known [20] only 
for the parent-independent mutation (PIM) model, in which the probability of mutating from 
allele j to allele i depends only on the child allele i. For a general non-PIM mutation model, 
finding an exact, closed-form sampling formula has remained a challenging open problem. 

In this paper, we make progress on this problem by deriving approximate, closed-form 
sampling formulas that are highly accurate when the mutation rate is low. More precisely, 
given a sample configuration n and the model parameters (mutation rate and transition matrix 
P), we consider the Taylor expansion of the sampling probability q{n \ 9, P) about 6* = 0. As 
discussed later, if P is irreducible when restricted to the observed alleles in the sample, then 
the leading order term in the expansion is proportional to ^^I'^il^i, where is the number 
of distinct observed alleles in the sample configuration n. Hence, 

q{n \e,P) = 6'l°"l-iO(n | P) + 0(6ll°"l), (1) 

where Q{n \ P) is the leading order coefficient that depends on the mutation transition 
matrix P but not on the mutation rate 9. In this paper, we consider the problem of obtaining 
exact closed-form formulas for Q{n | P). As many organisms typically have small per-base 
mutation rates, our results are of biological interest. 

By restricting the set of events in the coalescent genealogy for a given sample, Jenkins and 
Song [12] provided closed-form formulas for Q{n \ P) for an arbitrary transition matrix P 
when \On\ < 3. In this paper, we provide new proofs of those results, and extend them by 
supplying a closed-form formula for Q{n \ P) when = 4 and the transition matrix P is 
reversible restricted to the observed alleles. We prove our results using martingale arguments 
and use an urn construction related to the coalescent to develop a recursion for the approximate 
sampling probability, which can then be solved in closed-form using combinatorial techniques. 
As a corollary of our results, it can be seen that the simple general formula in [12, Theorem 
6.3] for Q{n \ P) when P is parent-independent restricted to the observed alleles also holds 
when P is reversible restricted to the observed alleles, provided that \On\ < 3. That formula 
fails to hold when |0„ | = 4 and P is not parent-independent restricted to the observed alleles. 

As there are four distinct DNA bases, our extension to the \On\ = 4 case seems natural. 
A more interesting reason is as follows: In multi-locus models with finite recombination 
rates, no closed-form sampling formula is known, even for the simplest case of two loci 
with either infinite-alleles or finite-alleles PIM models. However, recently a new framework 
based on asymptotic series has been developed [2, 10, 11, 13] to derive useful closed-form 
results when the recombination rate is moderate to large. The main idea behind that research 
is to perform an asymptotic expansion of the sampling probability in inverse powers of the 
recombination rate. We note that our one-locus sampling formula for the = 4 case 
provides an accurate approximation of the sampling probability for a completely linked (i.e., 
with zero recombination rate) pair of loci with two observed alleles at each locus (as is 
typical in single- nucleotide polymorphism data). Hence, our work serves as a starting point 
for finding approximate two-locus sampling formulas when the recombination rate is small, 
complementary to the earlier work [2, 10, 11, 13] for large recombination rates. We leave this 
problem for future research. 

We remark that, for a given sample configuration n and fixed parameters 9 and P, the exact 
sampling probability q(n \ 6, P) can be found numerically by solving a system of coupled 
hnear equations in 0(|n|^) variables, where \n\ denotes the total sample size and K denotes 
the number of allele types in the assumed model. One of the main motivations of our work is 
to remedy this high computational complexity. Evaluating our closed-form approximations is 
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much more efficient, in both time and space complexity. 

The rest of this paper is structured as follows. In Section 2, we lay out the model and no- 
tation used throughout the paper. In Section 3, we summarize our main closed-form sampling 
formulas, which we prove in Section 4 using martingale arguments and an urn construction. 
Numerical experiments demonstrating the usefulness of our approximate sampling formulas 
are provided in Section 5. 

2. Model and notation 

We consider Kingman's coalescent with a JC-allelic recurrent mutation model specified by 
the population- scaled mutation rate 6/2 and ergodic transition matrix P, where Pji denotes 
the probability of allele j mutating to allele i forward in time given that a mutation occurs. The 
stationary distribution of P is denoted by tt = (tti, . . . , ttk)- 

The following definitions will be used throughout: 

Definition 1. {n, sample configuration.) A sample of individuals is denoted by n = {ni)i^[K], 
where G Z>o denotes the number of individuals in the sample with allele i. The size \n\ 
of the sample n is denoted by the same letter in non-bold-face, n. For notational convenience, 
we use e, to denote the sample configuration with a single individual of type i and write 
n = n\e\ -\ h riKe-K- For a subset S C [K], we define ns = J2ies ''^i^i "■•5 ~ I'^sl- 

Definition 2. (On, observed allele types.) Given a sample n, let On C [K] denote the set of 
observed allele types; i.e., 0„ = {i £ [K] \ rii > 0}. The number of observed allele types is 
denoted by 

When the indices h,i,j,k and I are used in indefinite summations or products, they are 
assumed to range over On, unless stated otherwise. 

By exchangeability, the probabiUty of any ordered sample with configuration n is invariant 
under all permutations of the sampling order. We use q{n \ 9, P) to denote the stationary 
sampling probability of any particular ordered sample with configuration n. From the standard 
coalescent arguments [7, 8], it can be deduced that q{n \ 9,P) is the unique solution to the 
recursion 

n{n-l + 0)q{n \ 0,P) = ^ni{ni-l)q{n-ei \ 9,P) + 0J2 Pjiniq{n-ei + ej \ 0,P), 

i i,j 

(2) 

with boundary conditions 

q{ei \0,P)= TTi, foralH G [K]. 

If P is irreducible when restricted to the observed alleles On, then by unwinding recursion 
(2), it can be seen that — 1 is the smallest power of with a non- vanishing coefficient 
in the Taylor series expansion of q{n \ 6, P) about ^ = 0. Intuitively, for a sample with m 
distinct observed alleles, the coefficient of 0"^~^ in the Taylor expansion corresponds to the 
total probability of coalescent genealogies with the most parsimonious number (i.e., m — 1) 
of mutations. That P is irreducible when restricted to On is a sufficient (but not necessary) 
condition for the existence of such a parsimonious genealogy for sample n. 

Letting Q{n \ P) denote the coefficient of 0\'^^\~^ in the Taylor expansion, q{n \ 9, P) 
can be written as in (1). For simplicity, in what follows we simply write q{n) and Q{n) instead 
of q{n I 0, P) and Q{n \ P), respectively. 
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We now introduce some notation used throughout the paper. For a sample configuration n, 
we define the combinatorial quantity A(n) as 



A(n) 



(n-1)! • 



(3) 



For k e Z>o, the fcth falling factorial of x (denoted and the fcth rising factorial of x 

(denoted {x)kf) are defined as 



with {x)oi = (x) 



ot 



{x)ki = x{x-l)---{x-k + l), 
{x)kf = x{x + 1) • • • (a; + - 1), 

1. The fcth harmonic number H/. is defined as 

1 



Hk 



+ 



k' 



with Ho = 0. Given a sample configuration n = (ni , . . . , hk), a if -tuple m = (mi , . . . , m^) 
satisfying ^ m ^ n means < rrii < rii for all i E On and rrii = for all i ^ On, while 
^ m ^ n means < mj < for all i £ On and = for all i ^ O^. Also, 
^ m ^ n denotes < m, < for alH € [-ft"]. 

3. A summary of closed-form results for Q{n) 

In the case of \On\ = 1, it is easy to see that Q{n) = Hi for n = nei. In this paper, 
we derive closed-form expressions for the leading order coefficient Q{n) when \On\ < 3 and 
P is an arbitrary mutation transition matrix that is irreducible when restricted to the observed 
alleles On', and also when \On\ = 4, and P is irreducible and reversible when restricted to 
On (i.e., TTiPij = TTjPji for all i, j e On)- These closed-form results are summarized below. 

Theorem 1. For \On\ = 2 and P an arbitrary mutation transition matrix that is irreducible 
when restricted to On, Q{n) is given by 



Q(n)=A(n) ^ ^tt.-P,,. 



Theorem 2. For | On \ = 3 and P an arbitrary mutation transition matrix that is irreducible 

when restricted to On, Q{n) is given by 



Q(n)=A(n) ^ 

clistincti,j,fc£0„ 



n{nj + Uk 
niUjUk 



rik + 1)34. 

rijUk 



1) n{ni + Uk) 

[Hn — -ffjii-l) 



n{nj + nk)2i 



+ 2- 



n^UjUk 



n{nj + rife - 1) n{nj + nk)2i 



- 2 



niUjUk 



{rij +nk + 1)34. 



Approximate sampling formulas for general mutation models 



5 



Corollary 3. Suppose \Or. 



3 with sample configuration n = UaCa + ribei, + ricCc 



where a, b, c are distinct alleles in [K]. If the mutation transition matrix P is reversible and 
irreducible when restricted to the observed alleles On, Q{n) is given by 

— T^aPabPac H T^bPbaPbc + —T^cPcaPcb ■ 

n n n J 

Theorem 4. For \On\ — 4, if the mutation transition matrix P is reversible and irreducible 
when restricted to the observed alleles 0„, then Q{n) is given by 

Q{n) =K{n) ^ [TTiPijPikPin{n,i,j,k,l)+TTiPijPikPjiS{n,i,j,k,l)], 

distinct iJj/CjieCi 



where 

7{n,i,j,k,l) 



rii - 1 



2n n 



J "I 



2{ni 



■ rii 



1) {rii + rij + nk) 



ni 



2{nj + rife + ni) 



mini - 1) 



2njni 



+ 



{nk+ni){ni + nj -1) (ni + 71^)2;. 
2ninjni _^ 2ninjni 

(nj + rij + rife + 1)34. " "'"^ (n? + rij + 1)34. 



{Hn Hnf.+ni — l)i 



and 



S{n,i,j,k,l) = 



+ 



2njni 



ni + rij+nk-l (nj + + rife) 24,. 

rijUi 2njni 



+ 



[nk + ni){ni + Uj - 1) (rij + 71^)24.. 
2ninjni 2ninjni 



-{Hn — Hni-l) + 



{rii + rij + nfe + 1)34. ' {rii + nj + 1)34. 

4. Proofs of the main results 



{Hn - Hn^+m-l)- 



In this section, we construct an urn process to derive the closed-form formulas for Q{n) 
mentioned in the previous section. We use the urn process to decompose Q{n) into a sum- 
product of two vectors, one which depends only on the sample configuration n and the other 
which depends only on the mutation transition matrix P. Using this decomposition, we show 
that Q{n) corresponds to the probabihty of a certain event in the urn process. 

Throughout, we use R{n) to denote the following rescaled version of Q{n): 



R{n) 



Q{n) 
A(n)' 



(4) 



where A(n) is the combinatorial coefficient defined in (3). 
4.1. Description of the urn process 

Let n be the sample configuration of interest. We have an urn with n balls, n, of which 
have color i. We remove balls one at a time uniformly at random until there are no more balls 
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in the urn. However, whenever we "kill" a color (i.e., remove the last ball of that color), we 
add back a ball of a different color. We do this by picking another ball from the urn, copying 
it, and returning both copies to the urn. Note that when we kill the last color, we do not add 
any balls back, since there are no more colors to choose from. 

Suppose that when we kill color i, we add back a ball of color j. We then call j the parent 
of i, and call the last surviving color the wot. This generates a rooted tree whose vertices 
consist of the \On \ observed colors (alleles). 

Let T be any rooted tree on On. We denote the probability of generating T under the above 
process as P„(T). Let E{T) be the edge set of T, and let p{T) denote the root vertex of T. 
By convention, we draw edges as pointing away from the root, so the edge [j i) indicates 
that j is the parent of i. 

The main idea of this section is that to compute Q{n), it is enough to compute ^n{T) for 
each T. In particular, we prove the following theorem in Section 4.2: 

Theorem 5. Recall that for a transition matrix P that is irreducible when restricted to On, 
Q{n) denotes the first nonzero coefficient in the Taylor expansion (1) of q(n) about = 0. 
Given a rooted tree T described above, define fp{T) as 

fp{T)=^piT) n p,.. 

(j^i)eE{T) 

Then, the quantity R{n) = Q{n)/A{n) is given by 

R{n) = ^¥n{T)fp{T) = En[fp{T)], (5) 

T 

where the sum is taken over all rooted trees T with |0„| vertices bijectively labeled by On- 
That is, R{n) is the expectation of fp{T) under the above process. 

Note that we can view fp{T) as a probability as well. In particular, suppose we relabel the 
vertices of T as follows: we assign a new label from [K] to p{T) according to the stationary 
distribution w, and for each edge in T, we assign a new label to the child according to the new 
label of its parent and the transition matrix P. Then fp (T) is the probability that we assign 
the original labels to all the vertices, given that we drew T. That is, if Co„ is the event that we 
assign the original labels to all vertices, then 

fp{T) = P(Co„ I T) = 7r,(T) n P^^- 

This inmiediately leads to the following interpretation: 

R{n) = ^P(Co„ I T)¥n{T) = P„(CoJ. (6) 

T 

That is, R{n) is the unconditional probability that we correctly label all the alleles, if we use 
the urn process to generate a tree on the alleles and then use the tree to assign labels. 

4.2. An inductive proof of Theorem 5 

In this subsection, we provide an inductive proof of Theorem 5. In Section 4.3, we provide 
an alternative proof based on a modified coalescent process which provides a more intuitive 
explanation for why the urn process works. 
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Proof of Theorem 5. Recall the recursion in (2): 

i{n-l + 0)q{n) = ^nj(ni - l)q{n - Ci) + 9^ Pj^7Hq{n - ej + Sj). 



n 



Recall also that if P is irreducible when restricted to On, q{n) has leading order power 
Q\On\-i Taylor series. Hence we get the following recursion for Q{n): 

n{n-l)Q{n)= ^ ni{n^ - l)Q{n - ei) + ^ ^ PjiniQ{n - Si + ej). 

i:ni>l i:ni = l 

Plugging in Q{n) = A{n)R{n) and simplifying gives us the following recursion for R{n): 
n{n—l)R{n)= ni{n — l)R{n — ei) + ^ PjinjR{n — Si + Cj). (!) 

i:ni>l i-.rii — l j-.jy^i 

A simple induction over \On\ and n shows that this recursion has a unique solution given 
the boundary conditions R{ei). So if we can show (5) when \On\ = n — 1, and then show 
that J^T^i^On I T)Pn{T) satisfies the recursion (7), then we will be done. The base case 
is trivial: when On = {a}, there is only one possible tree, T = {a}, with Pn(r) = 1 and 
P(Co„ I T) = = lime^o = = A{n)R{n) = it!(n). 

To show J2t ^((^Or^ I T)¥n{T) satisfies (7), we start by giving recursions for Pn(T) and 
IP(Co„ I T). Let z(i) be the parent of i in T, and let L{T) be the set of leafs of T (where the 
root is not considered a leaf). Conditioning on the first event in the urn process gives us 

i:ni>l ieL(Ty.ni=l ^ ' 

Furthermore, if i € L{T), we have 

P(Co„ I T) = P,(,),,P(Co„\{.} I T \ {i}). (9) 
Using (8) and (9), and collecting terms, we arrive at 
n(n-l)^P(CoJT)P„(T) 

T 

T Vi:ni>l ieL{T):ni=l 

= ^ n,(n-l)^P„_e,(r)P(Co„ |T) 

i:ni>l T 

+ E E^^-^"^E^"-e.+e.(mco„\o}m, 

where the sum over T' is taken over all rooted trees with vertex set On \ {i}- Hence, 
P(Co„ I T)¥n (T) satisfies (7). □ 
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4.3. Connection to the coalescent 

In this subsection, we motivate our urn process by drawing a connection to the coalescent. 
We then use this connection with the coalescent to provide an alternate proof of Theorem 5. 

Let H be a history of mutation and coalescence events on n labeled individuals, and let 
g('H) be the probability of H. Then we have 



It turns out that only histories with exactly \On \ — 1 mutations contribute to the leading order 
term of q{n); this is the observation also utilized in [12J. Furthermore, each history of choices 
in our urn process corresponds with a genealogical history of |0„| — 1 mutations. This provides 
the basic intuition for why the urn sampling scheme works. 

We start by providing a modified coalescent that generates a history H that is consistent 
with n and has exactly \On\ — 1 mutations. We then show that this modified coalescent is 
equivalent to our urn sampling process. Finally, we prove Theorem 5 by relating the modified 
coalescent with Kingman's coalescent. 

Consider the following modified coalescent process on our sample: 

1. Select allele i with probability rrii/m, where m is our current configuration of alleles. 

2. If rui > 1, choose a random pair in allele i to coalesce (so m is replaced with m — ej). 

3 . If = 1 , have the last individual of allele i mutate to allele j with probability mj/{m— 
1) (so m is replaced with m — Ci + ej). 

4. Repeat steps 1 to 3 until all individuals have coalesced. 

It should be clear that the modified coalescent only generates histories with exactly |0„ | — 1 
mutations, since each mutation kills an allele permanently. 

If we take an unordered view of our sample, then the modified coalescent is equivalent to the 
urn process, for they have the same initial configuration and transition probabiUties between 
configurations. In particular, when m,i > 1 we move from m to m — with probability 
rrii/m, and when rUj = 1 we move from m to m — ej + Cj with probability m^/ (771)24,. 
generate trees on On by drawing an edge (j i) whenever we make a transition from m to 
m — Bi + Bj, i.e. whenever there is a mutation from i to j. 

We now give a proof of Theorem 5, using the modified coalescent in place of the urn 
process: 

Alternative proof of Theorem 5. Let 'H be a coalescent history with exactly M mutations. 
Running time backwards from the present, we suppose that the ith mutation was from allele 
Ui to allele v,, and that the most recent common ancestor has allele p. We further suppose that 
J, is the total number of lineages at the time of the ith mutation. Then we have that 



q{n) = 





(10) 



'H consistent with n 





-1 



2 



, and the 



{n-i+l)(n-i+e) 
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Now, observe that 

9W..i.^...(np...) „,(„_,),n«,,,,_„ - <■« 

Therefore, the Taylor series for q{l-L) has leading power O'^ , with coefficient Q{H). 

Hence by (10), the Taylor series for g(n) has leading power and its leading 

coefficient is given by the sum of all QCH) such that H is consistent with n and has \On \ — 1 
mutations. 

For such an H, let P„ {H) be the probability of generating H under our modified coalescent. 
Then we have that 

on—l 

To see this, note that if our current sample is m, the probabihty that the next event is a 
coalescence on allele i with > 1 is 

m,- 2 2 



m mi{mi — 1) m{mi — 1) ' 

and if mj = 1, the probability that the next event is a mutation from allele i to allele j (where 
j ^ i) is 

m{m — 1) 

Multiplying the probabilities of the mutation and coalescence events in H, and noting that the 
numerator of each mutation term cancels with the denominator of a future coalescence term, 
yields the equation (12). 

Combining (11) with (12) yields 



Q(7^) = A(n)7r^( [] P.,„,)Pr.(7^) 



Now let ^{H) be the resulting tree on On if we draw an edge {j — )• i) when allele i 
mutates to allele j. Then we have 



Q{n) = ^(^) 

H consistent with n 
H has I On I — 1 mutations 

= A(n)5^7r,(T)( n E 

T (j-^i)(zT ■H:.9('H)=T 

= A(n)^7r^(T)( n Pi)^r^{T) 



and hence 



T 

as needed. □ 
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4.4. A martingale property 

Here, we prove Theorem 1 and Corollary 3 by using martingales to compute P„(r) for 
On = {a, b}, and for On = {a, b, c} when P is reversible when restricted to On- We run 
time as follows: whenever we remove a ball in the um process, count this as one time step. If 
in doing so, we kill a color, count the adding of another ball as a separate time step. 

Let be the a-algebra generated by all sequences of choices up to time t. Let Xt be the 
proportion of balls that have color a at time t; so Xq = Ua/n. It is easy to check that {Xt} is 
a martingale with respect to {Tt}: Suppose that m is the remaining sample after time t — 1, 
and we are deleting a ball at time t. Then, 

m m — 1 ^ m m—1 m 

On the other hand, if we are adding a ball at time t, then 

f I X- 1 mavna + l ^ yr^nii rua ma „ 
m m + 1 ^ m m + 1 m 

So, {{Xt,J^t),t > 0} is a martingale. 

Proof of Theorem 1. Suppose C„ = {a, b}. Let T be the tree whose vertex set is On, with 
a being the root. Let r be the the first time we kill a color. Noting that r is a stopping time, we 
obtain 

P„(T)=E[P„(T|7;)] 

= E[I(Color a is the last remaining at time r)] 

= E[Xr]=E[Xo] = ^. 

n 

Therefore, by Theorem 5, 

Q{n) = A{n)('^7raPab + -^bPha). □ 

\ n n / 

Proof of Corollary 3. Suppose On = {a, b, c} and P is reversible when restricted to On- 
Note that P{Co„ \ T) does not depend on how T is rooted, for by reversibility we can move 
the root around by 

T^pPpk = T^kPkp, 'ik&On,k^ p- 

Therefore, we redefine Pn(r) to be the probability of drawing the undirected tree T. We still 
have R{n) = J2t ^i^o^ I T)Fn{T), but now the sum is taken over undirected T. Now let 
T be the tree on {a, 6, c} whose interior vertex is a. We draw T if and only if a is chosen 
as the parent of the first color that we kill. So, letting r be the first killing time and noting 
Xr =P„(T I JV), we have 

7? 

¥n{T) = E[¥n{T I Pr)] = = E[Xo] = ^. 

n 

Therefore, by Theorem 5, 

Q{n) = A{n)('^naPabPac + -TTbPbaPbc + -T^cPcaPcb) - □ 

\ n n n ) 
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4.5. A recursion for R(n) 

In this section, we derive a recursion for R{n) which will be useful for deriving closed- 
form formulas for Q{n) when \On\ = 3, 4. Given a sample configuration n and a subsample 
m, define the expression (^^j) as 



n \ 1 r / ri; 



n 



mj yrrii 



The following proposition provides a recursion relating R{n) to R{m) where | = |0„ | — 
1. 

Proposition 6. Suppose P is irreducible when restricted to On and let 6^^"-^~^Q{n) = 
6'l^"I^^A(n)i?(n) denote the leading order term in the Taylor expansion (1) of q{n) about 
6 = 0. Then, R(n) for |0„| > 1 satisfies the recursion 

R{n)- 2. 2. -py — ^^(^^^31) — ' 

mi=l 

with boundary conditions 

R{n) = TTa, 

for all sample configurations n = Haea, where a G [K]. 

Proof of Proposition 6. We can derive this recursion from the urn process as follows. Let 
Dij (m) be the event where the first killing replaces a ball of color i with a ball of color j, and 
where m is the (unordered) configuration immediately before this killing. Then for any event 
A, 

¥n{A) =Y. Pn(Ai(m))Pn(^ I Ai(m)) (13) 

mi=l 

where we use the fact that Vn{Dij{m)) = if rrii 7^ 1 or mj = for any j G On- 

We compute P„(£'ij (m)) when m )^ and rrii = 1. The probability that m is the 
remaining configuration after n — m draws is 

(n-w)! Uki^k)nk-mU _ (m) 



Uki^k-rnkV. (n)„_„4, (^) ' 

To see this, note that the first term is the number of ways we can make n — m draws that result 
in the configuration m, and the second term is the probabiUty of each such sequence of draws. 

When our current configuration is m with = 1, the probability that on the next draw we 
replace the last ball of color i with a ball of color j is rrij/ {171)21- Hence we get that 

P„(A.-(m)) = 



□ m{m-iy 



when my and m, = 1. 
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Plugging this into (13) yields 



(A)-\^ V I 



,{A\Di,{m)). 



(14) 



Now recall from (6) that R{n) — P„(Co„). That is, R{n) is the probability that we assign 
the original labels to all alleles, if we use the urn process to generate a tree on 0„ and then use 
the tree to assign new labels to the alleles. Note that 

P(Co„ I Aj(m)) = Pji^m-e,+e,{Cor,\{i}) = PjiR{m - Bi + Bj), 

since we need to use the urn process with sample m — + ej to correctly relabel On \ {i}, and 
then assign the correct label to {«} with probability Pji. Plugging this into (14) with A = Co„ 
yields the desired recursion, 



D/„^ _ p V- (m) mjR{m-Bi + Bj) 

mi=l 



□ 



In the next two subsections, we use the recursion in Proposition 6 to provide proofs of 
Theorem 2 and Theorem 4. 

4.6. Proof of Theorem 2(1 e>„ I = 3) 

For \On\ = 3, the following expression for R{n) can be derived using Proposition 6: 



_ \^ p sr a mjR{m-e, + Bj) 



^ ^ ( )m{m-l) ■f-^ m 

ljtkmdk,ljti 



E E ^ (") m2(m-l) 



mi = l 



mj(mj + l)7rjPj7 + ^ mjmknkPkj 



= E^.^E E 

Zjj^i m=3 0-<m:<Ti: 

mi=l, |rn.|— m 



□ m^(m-l) 



= E E E ■ 

i,j,k distinct m— 3 O^m^-n,: 

mi = l,|m|=m 



^ mj{mj -\-l)7TjPjk-\- ^ mjmkTTkPkj 

k:k^i,j k:k^i.j 

im) '^jPjiPjkmj{mj + 1) +TrkPkjPjimjmk 



'm?{m — 1) 



(15) 
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where in the second equality, the formula from Theorem 1 is used, noting that \Om,-ei+ej \ = 
2. If we define the quantities a{n, i, j, k) and j3{n, i, j, k) as 



n fn\ 

■^^rn^fm— 1) 

m=3 ^ ' Q~^m:<n: ^"iJ 

mi = l,\m\=m 



mj{mj + l), 



(16) 



and 



n ^ 



m=3 



(") 



then (15) can be rewritten as 



R{n)= ^ TrjPjiPjka{n,i,j,k) + ^ TTkPkjPji0{n,i,j,k). (17) 

distinct distinct 

Now consider a{n, i, j, k) defined by (16). We can remove the restriction in the inner sum that 
TOj = 1 by defining m' = m — ej, and so \m'\ = m—1. Also, since j ^ im (17), m'^ = nij. 
Making this change of variables from m to m' in the inner sum of (16), we get 



-^m (m- + 1) = ^""^^ n- 

/n\"''J\"''3^'-) fn\ (n-ni\ ■■"JV'J 

mi=l,|Tn|=m |m'|=m-l 



E 



E 



^ ™' 'm'Am'. + l). (18) 



Using identity (31) in Fact 5 of the Appendix, the summation over m' in (18) can be written 
as 



/n—niei\ 

O^m'^n—riiei: ^rn—lJ 
\'m'\=m—l 



= E (-1) 



|T| 



TG[L]: 



K)24.(to - l)2i ^ 2nj{m - 1) 



(n - rij - nT)24. n-ni-nr. 



\ ni-1 ) 
/ n— 
\m-l) 



(19) 



The only sets T satisfying the conditions in the summation in (19) are T = and T = {k}. 
Hence, substituting (18) and (19) in (16), we have 



a{n,i,j,k) = ^ 



^ 2r TT E rSirrijimj + 1) 



m=3 



1 



0-<m^n: 
m,; — l.|rM|— m 



Eto2(to-1) (") 
m=3 ^ ' \mJ 



f n—niei\ 

E V^"^^K + i) 



Im' \ =m—l 



E 



/n-nA ( (nj+Tn.\ 
«Vm-l/ I V m-1 / 



m=3 



i,/2i / ^ ^nAm — 1) 
''^^ (m- 1)21 + 2 



K' + nk)2i 



rij + Uk 
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rii 



TO=3 



E 



n 



\ m ) 

V m / 



\ TO / 



_ (rij + nk)2i TO + 1 rij + Uk m + 1 
Applying Facts 1 and 3 in the Appendix to (20) yields 

(nM ("^n ("to) , 2n,n, 1 



20) 



a{n,i,j,k) = ^ 



TO — 1 

n. 



■Hi 

n 



("to') ("to') inj+nk)2i TO + 1 

{nj)2i rij+Uk rij 



n { {rij + nk)2i 



rii + nk 



-2 



rijUk 



{rij + nk)2i 



Uj + nfc + 1 



riirij 



- 2- 



niUjnk 



n{nj + Tifc - 1) n(ri., + rifc) + Tife)2; 



+ 2 



(rij + nfc + 1)34. 



(21) 



Following a similar line of computation as above, we can find a closed-form expression for 
/3(n, i, j, A:) as follows: 



1 



1 

P{n,i,j,k) = V 

m—c 
n 

= E 

m=2 
n 

= E 



(") 



m\m-l) E 



m^ — l,|m|=m 

/ n—ni\ 
\7n-l) 



E 



m2(m-l) ' ^ f'^-^;) 

|r7^' I — rn— 1 



/ n—niei\ 
^ ™' ^TO^TO'fc 



\ TO-1 / 



rijUk 



TO=3 



(m-1) □ '(nj + nfc)24 



(m - 1)24 



y- ^jHk ("^+"'°) 2_ 

^in(ni + nfc)24 ^ V m+l 



n [nj + nfc)2^ i n 



Uj + rife 



- 2 



n,j + nj, + 1 



(Hn — Hm-i) — 1 



rijUk 



+ 2 



niJijUk 



n{nj + rife - 1) n(nj + 71^)24. (n^ + rifc + 1)34. 



(iJ„ — il„._ij[?2) 



where the second equaUty above is the same change of variables from m to m' = m — as 
in the a{n, i, j, k) term. The third equaUty follows from identity (32) in Fact 5, and the second 
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to last equality follows from Facts 1 and 3. Substituting (21) and (22) into (17), and using (4) 

gives 

Q{n) = A{n) ^ [T^jPj^Pjka{n,i,j,k)+TTkPkiPiil3{n,i,j,k)] 

i,j,k distinct 



A(n) J2 YjPjiPjk 



i,j,k distinct 



ntrij 



n{nj + rife - 1) n{ni + rife) n{nj + nk)2i. 
(rij + rife + 1)34. 



rijUk 



+ 2- 



n{nj + rife - 1) n{nj + nk)2i 



'^7 , , IN {-tin - J^m-i) 

{rij + rife + 1)34. 



Note that if P is reversible when restricted to the observed alleles 0„, then (23) simplifies to 
the expression given in Corollary 3. □ 

4.7. Proofof Theorem 4 (I On I = 4) 

Using Corollary 3, we first note the following alternate expression for R{n) when = 3 
and P is reversible restricted to the observed alleles: 



i,j,k distinct 



(24) 



Suppose \On\ = 4 and assume that P is reversible restricted to the observed alleles On. Then 
using Proposition 6, we obtain 

Ti(nn\ V" P (^) m,,i?(m - + eft) 



mi = l 



^) ruh \- rrii + Si^h PijPik 

2^ TT,;^^ 



^ -P^' (") _ 1) 

l,hjtl 0^m:<n: ^mJ ^ ' i,j,fc distinct, 

mi=l i,j^k^l 



m 



sr K p. p. p sr (^) + 



i,j,k,l distinct 



mi=l 



Vm/ 



+ ^iPiiPikPji E 7^ 

i,j,k,l distinct 



- (")m2(m-l)' 

mi=l 



where the second equality follows from using (24) since P is reversible when restricted to the 
alleles {i,j,k} C On ■ Similar to the proof in Section 4.6, if we define quantities ({n, i,j,k,l) 
and 6{n, k, I) as 



C{n,i,j,k,l) = 



7Tn?(m — 1) 

mi = l,|m|=m 



Y j^m,(m, + l), 
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and 



S{n,i,j,k,l) = 



E 



ni^lm — 1) 

m=4 ^ ' 0~i.m<n 



(") 
mi=\,\m\=m 



miirij, 



then, using (4) and (25), we obtain the following expression for Q{n) = A{n)R{n): 

C{n,i,j,k,l) 



Q{n) = A{n) ^ 

i,j,k,l distinct 



TTtPijPikPil- 



T^iPijPikPjl^in, k, I) 



.(26) 



By a very similar calculation to that in Section 4.6, using Facts 1 and 3, and identities (31) 
and (32) in Fact 5 of the Appendix, we obtain the following closed-form expressions for 

C{n,i,j,k,l) and 5{n,i,j,k,iy. 

C{n,i,j,k,l) 

_ni(ni+nj+nk ("-j;)2i rii 



ni 



{iii + rij + nk)2i rij + Uk + m 



+ 



n 



[Hn - Hni-l) - 1 



2ni{nj +nfc) 
(rij + Uj + nk)2i \ni + rij + rik + 1 

Uk + ni (ui + nj)2i ' (rh + nj)2i \ni + Uj + 1 

n 



rii + rife {ni)2i _^ ^muk 



Uj + ni (rij + nfe)24. (n, + nk)2i. \ni+nk + l 



{Hn — Hn^+m-i) — 1 

{Hn — Hn,+ni-l) — 1 



and 

5{n,i,j,k,l) 



n 



+ Tij + rife 



fiirii 



2ninj 



{rii + Uj + rife) 24. 

n 



(rii + Uj + nk)2i \ni + Uj + Uk + I 



-{Hn — Hni-l) — 1 



rij + Uj TliTlj 



2ninj 



1 — 7 — 1 — ^ 7 — , — \ — — I rri^"' ~ Hn^+n,-i) — 1 

nk + ni {rii + nj)2i {rii + nj)2i \ni + + 1 

Simplifying the expression for S{n,i,j,k,l), we get the expression stated in Theorem 4. 
Observing that ({n,i,j, k, I) is synnmetric in j and k, we see that for all i, j, k, and / distinct 
in On, 

C{n,i,j,k,l)+C{n,i,k,j,l) 



: j{n,i,j,k,l) +'y{n,i,k,j,l), 



where 7(n, i, j, k, I) is given by: 

rij - 1 



'y{n,i,j,k,l) 



n 



2njni 



2{ni + rij + rife - 1) (rij + rij + nk)2i 

ni{ni - 1) 2njni 

{nk + ni){ni + rij - 1) (n, + nj)2x. 



+ 



ni 



2{nj + rife + ni) 
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(n^ + Uj +nk + 1)3; K + nj + 1)3; 

Using the fact that T^iPijPikPu is also symmetric in j and k, we can then rewrite (26) as 

Q(n)=A(n) ^ [7riPijPikPin{n,i,j,k,l)+TTiPijPikPjiS{n,i,j,k,l)]. □ 

i,j,k,l distinct 

5. Empirical study of accuracy 

Here, we investigate the accuracy of approximating the sampling probability q{n) by using 
only the leading order term 6'l'-'"l~^Q(n). In this study, we solve the recursion (2) numerically 
to obtain the true sampling probability q{n) for moderate sample sizes. 

For a given sample n, define the approximate sampling probability, gapprox(^i). by 

aapprox(n) =0l°"l-ig(n). 

We can then define the relative error, Err(n), of the approximation Q'approx('^) from the true 
sampling probability g(n) as 

^ \q{n) - gapprox(n)| 
q{n) 

For a given sample size n, another natural measure of the approximation quality is the expected 
relative error under the distribution arising from the coalescent on samples of size n. Since 
q{n) is the probability of a particular ordered sample consistent with n, the probabihty p{n) 
of the unordered sample n, when samphng order is ignored, is given by 

pin) = " ^ 

\ni,...,nKj 

We can then define the expected relative error for a sample size n by AvgErr(n), given by 
AvgErr(n)= ^ p(n)Err(n)= i „ J - aapprox(n)|. 

■n:\n\=n ■n:\n\=n 

We also define the worst-case relative error, WorstErr(n), for a given sample size n as the 
worse relative error among aU samples of size n. Specifically, 

w /- \ T3 ^ \ - 9approx(n)| 

WorstErr(n) = max Err(Ti) = max ^-^ -r-^ — ^ — -. 

n:\n\=n n:\n\=n Q{n) 

To study the accuracy of approximating q{n) by gapprox(»i)> we examine the behavior 
of AvgErr(n) and WorstErr(n) for a transition matrix estimated from real biological data. 
Specifically, we use the reversible phylogenetic mutation rate matrix estimated in [21, Table 
1, matrix (1)] for the V??-globin pseudogenes of six primate species. Since their estimated 
matrix is a matrix of nucleotide substitution rates used for phylogenetic analysis, we rescale 
it by the minimum amount that can make it a valid Markov transition matrix. This rescaled 
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150 200 250 
Sample size n 

(a) 




150 200 260 
Sample size n 

(b) 



Figure 1: Error plots as a function of the sample size n, for the transition matrix P in (27) 
and mutation rate 6 € {10~^, 5 x 10~^, 10~^}. (a) The expected relative error, AvgErr(n). 
(b) The worst-case relative error, WorstErr(n). 

matrix, denoted by P, is given below to three digits of precision, and is used in our numerical 
experiments with different values of the mutation parameter 6: 

( 0.433 0.398 0.074 0.095 \ 

0.665 0.000 0.164 0.171 

0.074 0.098 0.394 0.434 ' ^ ' 

\ 0.147 0.159 0.674 0.020 / 

in the (T, C, A, G) basis. The stationary distribution corresponding to this transition matrix is 
TT = (0.308, 0.185, 0.308, 0.199) to three digits of precision. 

For many neutral regions of the human genome, typical mutation rates per base are in the 
range lO'^ < 6 < lO'^ [16], and we consider 9 € {10-3,5 x 10"^ IQ-^} in our study. 
For the transition matrix in (27), the expected relative error AvgErr(n) and the worst-case 
relative error WorstErr(n) are plotted in Figures 1(a) and 1(b), respectively, as functions of the 
sample size n. As can be seen from the plots, both the expected relative error and the worst- 
case relative error grow very slowly with the sample size n. Further, the ratio of WorstErr(n) 
to AvgErr(n) is a small number between 1.3 and 2.1 for all n < 360, and is decreasing in 
n. Hence, it appears that the approximation quahty of gapprox(^i) is uniformly good over all 
samples n for any given size n. 
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Appendix 

Here, we provide some general combinatorial identities which are used several times for 
proving the main results in this paper. 
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Fact 1. For any positive integers x, y, a and b where b < a and x <y. 



y /6\ (a+l-x\ _ ( a-y \ 

E \m) _ \a+l-b) \a+l-b) 

=x Vm/ \6/ 



(28) 



m=x \m 

Proof. Starting from the left hand side of (28), we have: 



y 



E \mJ 

m=x \ml 



b\{a-b)\ 



al 



y 

E 



a — m 
a — b 



/a+l-x\ _ ( a-y \ 
_ \a+l-b) \a+l-b) 

where the last equahty follows from the standard combinatorial identity that for all positive 
integers a, n, and k. 



E 



n — I 
k 



n — a + 1 
k + 1 



□ 



Fact 2. For positive integers a and b, 

1 fa — m 



E 



Fact 2 can be verified by induction [4] or by the method of Wilf-Zeilberger pairs [5]. 
Fact 3. For positive integers a and b where b < a, 



(a) 



a+1 



(Ha 



+ 1 



Proof. Starting from the left hand side of (29), we have: 



H. 



a—b. 



f C) 1 



6!(a-6)! 
a! 

^ 6+1 

^ E 

m=2 

b+1 



E 

m=l 



- 1. 



1 



(29) 



a — b I m + 1 



a+1 — m\ 1 



a—b 



m 



1 

(!) 

1 

(!) 

a + 1 
6 + 1 



E 



a + 1 — m\ 1 



rn=l 
a + 1 

6+1 



m 



a — 6 

{Ha+l — Ha-b) 
{Ha+l — Ha-b) — Ij 



where the fourth equahty follows from using Fact 2. 



□ 



We also list some facts about the moments of a hypergeometric distribution which are 
appealed to several times in the paper 
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Fact 4. If a multivariate hypergeometric distribution is parameterized byn ~ (ni, n2, . . . , n^), 
where n = \n\, and a sample of size m, m = {mi,m2, . . . ,mL), is drawn from it, then for 
any t= (ti, t2, • • • , ii,) where ti>0 for all i,t= \t\ and t < n, 



E 



til 



.i=l 



|m| = 



(30) 



Proof. Starting from the middle term in (30), we get: 



\m) 1=1 



\m.\—m 



E nf^l('^»)ti4- im-t) 
Y:;^ {m)ti -j—pY 



Iml — m 



\m-tJ 



2^ -7^ 

^^'^^ 0-<m-<n-t: I » 



\m\=m—t 



(n) 



4 



where the last equaUty follows because the term being summed is the probabihty mass function 
of a multivariate hypergeometric distribution parameterized by n — t, and the summation is 
over the entire domain of the distribution, and hence is 1 . □ 

In the following fact, we compute some second moments of the hypergeometric distribution 
parameterized by n when restricted to those samples m which are non-zero at aU types. 

Fact 5. If n = (ni,n2, . . . where n = \n\, and 1 < j ^ k < L, then we have the 
following identities: 



Iml—m 



TC[L]: 



{n-nT)2i n-riT 



tn-nT\ 
\ m ) 

(") 



^mjm,= E (-1)' \n-n,)2 n 

|m|=m j^T 



(31) 



(32) 



Proof. Applying the inclusion-exclusion principle and using Fact 4, the identity in (3 1) can 
be obtained as 



E S|-.K + i)= E (-1)'"' 

^•m' TC[L]: 



/n-nT\ 
\ m J 



\m\—rri 



/n-nT\ 
Km/ 

(") 



E (-1)'"' 



TC[L]: 



{nj)2i{m)2i _^ 2njm 



{n-nT)2i n-riT, 



\ m ' 

(") 

Km/ 
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Similarly for (32), we have 

Km/ , 



E 



0-<Tn:<n: 
|m|=m 



(") 



mjimk 



|T| 



TC[L]: 
3MT 



n—nT\ 
m I 



^ y (n—riT 
O^m^n—TiT' 

\m\—m 



fn-nT\ 
Km) 



m^rrik 



\ m ) 
Vm/ 



^ (n-nT)2i (") ■ 



□ 
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